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Abstract 
> ■ 

' We consider a reformulation of QED in which covariant Green functions are used to solve 

I for the electromagnetic field in terms of the fermion fields. It is shown that exact few-fermion 

' eigenstates of the resulting Hamiltonian can be obtained in the canonical equal-time formalism 

I for the case where there are no free photons. These eigenstates lead to two- and three-body 

Dirac-like equations with electromagnetic interactions. Perturbative and some numerical solu- 
tions of the two-body equations are presented for positronium and muonium-like systems, for 
various strengths of the coupling. 



o 
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' 1 Introduction 

It has been pointed out in previous publications ^ that various models in Quantum Field Theory 
(QFT), including QED, can be reformulated, using mediating-field Green functions, in such a way 
^ . that exact few-particle eigenstates of the resulting partially truncated Hamiltonian can be obtained. 

^ I This approach was then applied to two-body eigenstates in the scalar Yukawa (Wick-Cutkosky) 

theory |^, We implement such an approach to QED in this paper. 

The Lagrangian of two fermion fields, tp{x) and interacting electromagnetically, is 

C = i^{x){i-r^d^-qi-r^A^{x)-mi)'il;{x)+'${x){i'y^'d^-q2-f^A^{x)-m2)cl){x) 

- ^ {d^Afsix) - dpA^{x)) (d^A^ix) - d'^A'^ix)). (1-1) 

The corresponding Euler-Lagrange equations of motion are the coupled Dirac-Maxwell equations, 

(iT'^a^ -mi)V(x) = qi 7M^(x)^(x), (1-2) 
( i 7^9^ - m2)0(x) = q2 7''^m(x)(/)(x), (1-3) 

and 

d^d^A'^ix) - d'^d^A^'ix) = fix), (1-4) 
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where 

f{x) = qi tp{x)Ytp{x) + q2 (p{x)-f''(j){x). (1-5) 

The equations ( [L-2| ) - ( |1-4|) ca n be decoupled in part by using the well-known 0, ^ formal solution 
of the Maxwell equation (|1-4| ) , namely 



A^{x) = Alix) + / D,,,{x - x') f{x')d*x' 



:i-6) 



where D^^(x — x') is a Green function (or photon propagator in QFT terminology), defined by 



dad"D^^{x - x') - d^d"Da„{x - x') = g^u5^{x - x'), 



(1-7) 



and AP^{x) is a solution of the homogeneous (or "free field") equation {1^^ with j^{x) = 0. 

We recall, in passing, that equation (1-7) does not define the covariant Green function 
x') uniquely. For one thing, one can always add a solution of the homogeneous equation (Eq. ( |1-7|) 
with g^jj — > 0). This allows for a certain freedom in the choice of D^j^^ as is discussed in standard 
texts (e.g. refs. ^). In practice, the solution of Eq. (1£7), like that of Eq. (1^), requires a 
choice of gauge. However, we do not need to specify one at this stage. 

Substitution of the formal solution (l^O) into equations (1^2) and (1^3) yields the "partially 
reduced" equations. 



and 



{\rd^. - mMx) = gi7^ (^Alix) + J d^x'D,,,{x - x')f{x')^ ^{x), 
{ird^, - m2)m = q2l^ {aI{x) + j <fx'D^,{x - x')f{x')^ (t>{x). 



:i- 



:i-9) 



These are nonlinear coupled Dirac equations for two different fermion fields. To our knowledge 
no exact (analytic or numeric) solutions of equations (1-8) and ( |1-9| ) for classical fields have been 
reported in the literature, even for the case of a single fermion field (say = 0), though approximate 
(perturbative) solutions have been discussed by various authors, particularly Barut and his co- 
workers (see refs. ||6|, 0] and citations therein). However, our interest here is in the quantized field 
theory. 

The partially reduced equations ( |l-8|) -(|l-S|) are derivable from the stationary action principle 



6 Slip, (l)]=5 



J Cr (fx = 0, 



(1-10) 



with the Lagrangian density 

£.R = %|){x){\-f^^^-ml-ql'y^A'^{x))i){x) + 4){x){\-f^^^-m2-q2l^lA'^{x))(t){x) 

-\ J d'x'f{x')D,,{x - x')f{x), (1-11) 

provided that the Green function is symmetric in the sense that 

Dfj^y{x - x') = D^i,{x' - x) and D^y{x - x') = Du^{x - x'). (1-12) 

One can proceed to do conventional covariant perturbation theory using the reformulated QED 
Lagrangian ( 1-11 ). The interaction part of ( 1-11 ) has a somewhat modified structure from that 



2 



of the usual formulation of QED. Thus, there are two interaction terms. The last term of (|T^ 
I 11| ) is a "current-current" interaction which contains the photon propagator sandwiched between 
the fermionic currents. As such, it corresponds to Feynman diagrams without external photon 
lines. The terms containing corresponds to diagrams that cannot be generated by the term 
containing D^^, particularly diagrams involving external photon lines (care would have to be taken 
not to double count physical effects). However, we shall not pursue perturbation theory in this 
work. Rather, we shall consider an approach that allows one to write down some unorthodox but 
exact eigenstates of a truncated model, in which terms involving Aq are ignored. 

The paper is organized as follows. In section 2 we quantize the system using the canonical equal 
time formalism in the Schrodinger picture. In section 3 an unconventional "empty" vacuum state 
is used to construct exact one-, two-, and three-fermion eigenstates of the Hamiltonian, truncated 
to exclude states with free (physical) photons. In section 4 we show that the resulting two-fermion 
equation is the Breit equation in the Coulomb gauge, but that it is the Eddington-Gaunt equation 
in the Lorentz gauge. In section 5 we demonstrate that the Breit equation can be obtained in the 
Lorentz gauge, provided that higher-order retardation effects are taken into account. 

The reduction of the Breit equation to radial form is described briefly in section 6. For states 
of zero total angular momentum (J = 0), four coupled radial equations are shown to arise. The 
analytical structure of their solutions is studied in section 7. Perturbative O(a^) corrections to the 
Rydberg spectrum of J = states are obtained in section 8. In the case of equal rest masses the 
J = 0+ state equations have no unusual singularities and can be solved numerically. Some of these 
results are presented and discussed in section 9. The remainder of the paper is devoted to the study 
of J > states. In the section 10 the set of eight coupled radial first-order differential equations is 
reduced to four first-order ones and then to two second-order Schrodinger- like equations. They are 
solved perturbatively in section 11 and O(a^) relativistic corrections to the non-relativistic mass 
spectrum are obtained. A summary and concluding remarks are given in section 12. 

2 Hamiltonian in the canonical, equal-time formalism 

We consider this theory in the canonical, equal-time formalism. To this end we write down the 
Hamiltonian density corresponding to the Lagrangian ( p.-llj ) 

Hr = ij^x){-ia-V + miP)'ilj{x) + qi%l^{x)'jf,A'^{x))'il^{x) 

+ 0t(x)(-i«- V + m2/3)0(x) + g2^(x)7^A[^(x))</)(x) (2-1) 
+ ^Jd'x' j^^{x')D^,{x-x')f{x) 

where we have not written out the Hamiltonian density for the free Aq{x) field. 

Equal-time quantization corresponds to the imposition of anticommutation rules for the fermion 
fields, namely 

{iPa{x,t),i;l{y,t)} = {(l)aix,t),4'l{y,t)} = dap5'\x - y), (2-2) 

and all others vanish. In addition, if Aq ^ 0, there are the usual commutation rules for the Aq 
field, and commutation of the Aq field operators with the and (j) field operators. 

The Hamiltonian ( |2-1D contains an interaction term that is nonlocal in time, which can compli- 
cate the transition to a quantized theory. We shall avoid this problem by working in the Schrodinger 
picture with t = in the expressions for the field operators and currents, that is ip{x) = ip{x, t = 0), 
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j^{x) = j^{x, t = 0), etc. in equation ( |2-1| ). This corresponds to neglecting higher order retardation 
effects. Thereupon we obtain the result 

J dt' D^^{x - x') = G^^{x - x'), (2-3) 

where 

G^Ax) = J ^G^.(/e)e'^-^ and G^,{k) = D^^ik^^ = (0, fc)). (2-4) 
For example, in the Lorentz gauge {d^A^ = 0), we have 

G^,u{x) = g^,u—^. (2-5) 
47r kc 



Thus, in the Schrodinger picture, the third term of the Hamiltonian density ( |2-1[ ) takes on the form 

ni{x) = \j d^x' j>'{x')G,,{x - x')f{x). (2-6) 

In the remainder of this paper, we shall consider the simplified model without the interaction 
terms in (2-1) that contains A^. Such a model is suitable for describing few-fermion states in- 
teracting via virtual photon exchange, but without decay or annihilation involving free (physical) 
photons. In short, in all that follows we consider the field theory based on the Hamiltonian den- 
sity of ( p-l| ) but with Aq{x) = 0. An attractive feature of this model is that exact few-fermion 
eigenstates of the Hamiltonian can be obtained. 



3 One, two and three fermion eigenstates 

We consider now the model for which the Hamiltonian, in the Schrodinger picture with t = 0, is 
given by the expression 

Hr = Hjjj + H(j, + Hj, (3-1) 

where 

= J d^x (a;, 0)(— i Q! • V + mi/3)'0(a;, 0), (3-2) 

H^ = J d^x (l)'^{x,0){-ia-V + m2l3)(l){x,0), (3-3) 

and Hf = J d^xl-Li{x), where 'Hi{x) is given in equation ( |2-6D . Note, again, that the terms in Aq 
have been suppressed, so that processes in which free (physical) photons are emitted or absorbed 
are not accommodated. 

The Hamiltonian Hn has the same structure as the Coulomb-QED (CQED) Hamiltonian, that 
is the Hamiltonian of QED in the Coulomb gauge, but with the "transverse-photon" part (that 
contains a ■ A) turned off. Indeed Hj^ would beidentical to Hcqed if the indeces /i and v took on 
only the value in equation ( p-6| ) (as it is, fi,^ = 0, 1, 2, 3 in Eq.( |2-6| )). It has been shown earlier 

that exact two-fermion eigenstates of Hcqed can be written down if we use an unconventional 
(or "empty") vacuum state, |0), defined by 

ijaix,0)\0) = Mx,0)\6) = 0. (3-4) 
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The same is true of the present more reahstic model, as we point out below. 

The unconventional empty vacuum definition ( ^-4] ) means that Tp{x) is interpreted as a (free) 
Dirac-particle annihilation operator, while ^jJ^ (x) is, correspondingly, a Dirac-particle creation op- 
erator. By "Dirac-particle" we mean one described by the full Dirac spinor, including positive 
and negative frequency components. (Recall that in the conventional approach, i.e. using a Dirac 
"filled negative energy sea" vacuum which is annihilated by the positive frequency component of 
ijj, it is only the negative-frequency component of tp that is an antiparticle creation operator, and 
the positive-frequency component of that is the particle creation operator). 

With these conventions, we write the normal-ordered Hamiltonian 



iHr; = + H^ + ^ J d^x (fx' G^^{x - x') 



(3-5) 



where = ^(a^) and cj) = 4>{x'), etc. The normal ordering is achieved by using the anticommutation 
rules ( p-2| ) as usual; but note that it is not identical to the conventional normal ordering because of 
the unconventional empty vacuum that is being used, and the unconventional definition of ^l' and 
(j) as annihilation operators and of as creation operators. To underscore this unconventional 

procedure we use the notation ; Hr; rather than : Hji : in equation ( p-5| ). 
We note that the state defined by 



|1) = J d-^x'il^'{x)F{x)\0), (3-6) 

where F{x) is a 4 x 1 c-number coefficient vector, is an eigenstate of ;H^; (Eq. ( |3-5| )) provided 
that F(x) satisfies the equation 

{-ia-V + miP)F{x) = EF{x), (3-7) 

which is the usual time-independent one-particle Dirac equation (with positive and negative energy 
solutions), so that F{x) is a Dirac spinor. Therefore, we refer to |1) as a one-Dirac-fermion state. 
Similarly, the two-Dirac-fermion state, 

|2) = J d\d^yF^p{x,y)ijl{x)<t>l{ym, (3-8) 



(summation on a, /3 = 1, 2, 3, 4 implied) is an eigenstate of ; Hr; (equation (|3-5|)) provided that the 
4x4 eigenmatrix F satisfies the equation 

hmAx)F{x,y) + [hmMF^ix,y)Y + qiq2G^u{x - y)TF{x,y){fY = EF{x,y), (3-9) 

where hmj{x) = — ia • + mj/S, 7^ = 7o 7^ = (IjCt); ^-nd the superscript T indicates the 
transpose of the matrix in question. 

The detailed form of the interaction matrix G^y{x — x') depends on the choice of gauge. This 
equation (|3-S|) is a two-fermion Dirac-like, or Breit-like, equation with positive and negative energy 



solutions, and is, in this respect, different from those obtained in the conventional approach |9|-|15] 
in which the negative-energy solutions do not arise. 
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We note, in passing, that if the interaction is turned off in Eq. ( ^-9|) (i.e. (?i = ^2 = 0), then 
the solution can be written as 

F{x,y) = f{x)g'^{y), (3-10) 



where f{x) and g{y) are solutions of the one-body Dirac eigenvalue equation (|3-7| ). This indicates 
that m. F = [Fij], the index i corresponds to particle 1 (with coordinates x) while j corresponds 
to particle 2 (with coordinates y). 

In the rest frame of the two-fermion system (i.e. when |2) is taken to be an eigenstate of 
the momentum operator for this QFT, with eigenvalue 0), equation ( ^-9] ) reduces to an analogous 
equation in the single relative co-ordinate r = x — y: 

hrnAr)F{r) + [hm,{-r)F^{r)Y + qiq2G^,{r)rF{r){rf = EF{r). (3-11) 

It can, therefore, be reduced to a set of ordinary, coupled, first-order differential equations for states 
of given . Such equations can, at the very least, be solved numerically. This is a straightforward. 



though somewhat tedious, problem [Q, 16| which we address below. 

The structure of the Hamiltonian (Eq. ( |3-5D ) is such that generalizations to systems of 

more than two fermions are readily obtained. Thus, the three fermion state, corresponding to a 
system like |e~e^/i+), defined by 

(fxiSx2d?x:i Fa^a2a3{xi,X2,x:i) i;\^{xi)^l)\^{x2)(t)l^{x:i)\^) , (3-12) 

is an exact eigenstate of ; Hr; with eigenvalue provided that the 4^ = 64 coefficient functions 
-^ai«2a3(^ii a^s) satisfy the three-body Dirac-like equation, 

mi (332)]q,2q. acta (331 , 3^2, 333) 

+ [hm2{xz)\aia Fa^a2a{xi,X2,X^) + qiq2G^y{x-i - Xi) i^^)aia iY)asl3 Faazpixi, X2, X3) 

+ qiq2Gf,u{x3 - X2) {T)a2a i^") X2,X3) (3-13) 

{Xi -X2) i^^)aia {Y)a2P Fal3a3{xi,X2,X3) 
— F Fo^-^Q^^a^ {Xli X2 ) 333) , 

where summation on repeated spinor indices is implied. In the rest frame of the three-body system 
equation ( 3-13| ) reduces to one in two independent vectors only. Nevertheless, the reduction of the 



equation for states of given is more formidable than in the two-body case. Even then one is left 
with the full complexity of a relativistic three-body system. We shall not consider solutions of the 



three-body equation (3-13) in this paper 



4 Two-body equation in the Coulomb gauge 

At this stage we must specify the choice of Green's function, that is a choice of gauge. We could 
use any gauge, in principle, but we shall use the Coulomb gauge, as this avoids nonphysical degrees 
of freedom and the need to take account of auxiliary conditions (such as d^A^ = in the Lorentz 
gauge). We use the relation 
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and note that 



Dooik) = Doj{k) = 0, A, (A:) = ^[Sij 



Ivikj 



in the Coulomb gauge [1^]. Therefore, if we use Eq. ( |2-3| ) and the identity ([18|, Eq. 39. 



d?k : I, ^ a ■ k b ■ k 

e 



1 1 / a ■ r b ■ r 
a b 7. 



(27r)3 " (A;2)2 2r 

we obtain the coordinate-space representation of the Coulomb-gauge Green function 



Goo(r-) = — , Go. = 



1 



X 2 



(4-2) 



(4-3) 



(4-4) 



where r = |r| and r = (xi,X2,X3). Consequently, in the Coulomb gauge, equation ( |3-11 ) becomes 



hmAr)F{r)+ h^,{-r)F^{r) ' +V{ 



T 



F{r) 



cxF{r) ■ OL^ H — ■ cxF(r) r ■ a 



T 



47r r 



EFir), 
(4-5) 



where V{r) 

Keeping in mind the notation, as explained below Eq. ( 3-10| ), we see that equation ( [4-5D is 
nothing other than the Breit equation written in the rest frame of the two-fermion system. 



If we had used the Lorentz- gauge form of the Green function ( |2-5| ), equation (3-11) would take 
on the form 



h^,{r)F{r) + h^,{-r)F^{r) + V{r) F{r) - aF{r) ■ = EF{r), 



(4-6) 



which is recognized to be the Eddington-Gaunt equation |20, The Gaunt equation, unlike the 
Breit equation, does not contain even lowest-order retardation effects (see, also, section 5, following). 
Therefore, it will yield energy eigenvalues that will differ from those of the Breit equation already 
at O(a^). 

Note, also, that equations ( ^4-5| ) and ( [4-6| ) , with the terms involving a left out, become identical 
to the Coulomb-QED model discussed earlier H. 



5 Two-body equation in the Lorentz gauge 

Although we shall use the Coulomb gauge in this paper, it is instructive to see how, in the Lorentz 
gauge, one needs to keep retardation effects at least to lowest non-vanishing order in order to 
achieve the same results. The Green function of the d'Alembert equation in the Lorentz gauge has 
the form: 

D^^uix) = g^uD{x), D{x) = ^K^^)- (5-1) 
The reduced Lagrangian Cr in this case reads (with Aq{x) = 0): 

i^R = + + Ci 

= ii{x)(\-i^'^^,-ml)^l){x)+l^[x)(■l-f^'^^-m2)cp(x)-]^A^'{x)3^{x), (5-2) 
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where the potential of electromagnetic interparticle interaction is 

A^{x) = J <fx'D{x - x')f{x') = J Sx' j dt'D[{t - t'f -{x- x'f]f{t', x'). (5-3) 

Thus, the Lagrangian ( ^-2| ) is nonlocal in time. Because of this the standard hamiltonization 
procedure is not applicable. 

In order to employ the canonical Hamiltonian formalism it is necessary to convert this La- 
grangian to single-time form. We shall do so by employing a procedure which takes into account 
the retardation effects approximately. 

Using the substitution t' = t + X in (5-3) and expanding the current j in a Taylor series in A, 
we obtain the result 

jix',t') = j{x', t + X)= j{x', t) + Xjix', t) + ^X^j{x',t) + ■■■, (5-4) 

which reduces the potential ( |5-3| ) to the form: 

A^{x,t) = J d^x'!^G{r)f{x',t) + ^GMf{x',t) + ---Y (5-5) 

where r = \r\ = \x — x'\. The functions 

G{r) = [dX D{X^ -r^) = — (5-6) 
J 47rr 

and 

Gi(r) = / dXDiX^ - r2)A2 = ^ (5-7) 
J 47r 



G[{r)=rG{r). (5-8) 



satisfy the relation 

The terms of odd power in A in the expansion ( |5-4| ) vanish, because D is an even function of A (c.f. 
(^-6|), (|5-7|)). As a result, the interaction Lagrangian Cj (up to surface terms) takes the following 
single-time form: 

Ci ^ Cf^ + C^j^ + Cf\ (5-9) 



where 



= -\j d^x'G{r)f{x)f{x'), (5-10) 
4'^ = \j d^x'G{r)3{x).j{x'), (5-11) 



Cf = \j d^x'GMf{x)'f{x') (5-12) 

(hereafter we omit the common time argument t). 

The quantized theory based on the Lagrangian £(0) = + + cf^ (the Coulomb QED 
mentioned in section 3) was discussed in |^]. It takes into account the relativistic kinematics of 
the fermion fields exactly, but describes their electromagnetic interaction with the "transverse- 
photon" part turned off. The terms C^p and cf'^ can be treated as first-order corrections to £^^\ 
thus providing the approximate single-time form £5 = C'^^^ + + cf'^ for the reduced nonlocal 
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Lagrangian {c.f. Eq. ( |l-llD ). Other terms following from the expansion are corrections of 
higher order. They will not be considered in the present paper. 

The Lagrangian Cs leads to Euler-Lagrange equations which are second-order in time deriva- 

(2) 

tives, because of the term £} . Thus, it describes the system with twice as many degrees of freedom 
as jC^*^-* does, because ijj^ are no longer the conjugate momenta of ifj. This changes completely the 
dynamical content of the fields ip and (j). Since the second-order time derivatives occur in small 
terms only, they should be eliminated by means of the Euler-Lagrange equations of a lower-order 
approximation. But the resulting field equations are then not necessarily the Euler-Lagrange equa- 
tions of a known Lagrangian. Thus the transition to the Hamiltonian and to a canonical quantum 
description becomes unclear. 

To avoid this difficulty it is tempting to eliminate the time derivative of the charge density 

(2) 

directly in C\ by taking into account the conservation law 

d^f = Q, i.e. jo + V-j = 0. (5-13) 

This conservation law is a consequence of the Euler-Lagrange equations, which follow from the 
reduced Lagrangian Cr as well as from the truncated one C^^\ However, the direct use of the 
equations of motion (or their consequences) in the Lagrangian is not a correct procedure: it 
changes the equations of motion themselves. This fact was first emphasized in the case of the 



Golubenkov-Smorodinskii Lagrangian |22, 23] and then subsequently discussed in the literature 
[p^]-|28[|. Instead, one can use the method of "double zero", used in refs. p^ , p7[ . In our case this 
consists of the following modification of the Lagrangian: 

Cs^Cs = Cs + cf\ (5-14) 

where 

4'^ = -\j d?x'Gi{T) {j\x) + V • 3{x)] [f{x') + V • 3{x')] . (5-15) 
It is easy to see that the term Cj possesses the property: 



= 0, (5-16) 

<5 / d3xC(0)=Q 



so that it does not change the variational problem to the accuracy required. On the other hand, it 
cancels those terms of Cs which are quadratic in time derivatives of the fields. Thus the modified 
Lagrangian Cs yields equations of motion, which are first order in the time derivatives of the 
particle fields ip, (p. 

Next, we perform the following transformation of the field variables: 

= (1 - iq^W)i^ ^ e-'^^^Tp, V' = (1 + iqiW)'^ ^ e^'^'^i^, 

0^0= (1- ig2iy)0«e-*^2W^^ ^^^= (1 + i^2W^)^~ e^^^w/^^ (5.17) 

where 

Wix) = ^J d^x'Gi{r)V' ■ j{x'). (5-18) 

The transformation ( ^-ITp can be regarded as an approximate U{1) gauge transformation, which 
however is not canonical due to the dependence of W on the fields. This transformation removes 
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time derivatives from the interaction part of £5. To the accuracy required, and up to surface terms, 
Lagrangian £5 can be written in the form 



where the notations C^'^\ j, etc. mean that the fields ip, cp are replaced by ip, (p. 

Integrating the last term of £5 by parts, omitting surface terms, and using equations 
(^-8|) , we reduce the Lagrangian £5 to the form 



£5= i(^\^ + cP^^)-Us, 



where 



Us 



^'f'(a;)(-iQ! • V + miP)'ip{x) + ^ {x){-\cx ■ V + 'm2P)^{x) 
+ d'x'l [fix) fix') - \iix).3\x') - l-^{r . 3\x))ir . jix'))^ . 



(5-19) 



and 



(5-20) 



(5-21) 



This formulation allows us to treat the variables and cp^as the canonical conjugates of V' and 
(p, respectively. That is, we impose the anticommutation relations ( ^-2| ) on the underscored fields, 
and not on the original ones, when performing the quantization. Thereafter, since the Hamiltonian 
H Q = J cfixTig is formally identical with the Coulomb gauge Hamiltonian (c.f. (2-1) with (|4-4D ), 



calculations like those of section 4 lead to the Breit equation (4-5) 



6 Two-body equations in block component form 

Although two fermion equations have been around since the 1920s, their full reduction to radial 
form is of more recent vintage (see, for example |pg| ], |^]). The reduction of equation (|4-5D to radial 
form is essentially the same as presented in ref. [S], hence all the details shall not be repeated here. 



As shown in |8[, we note that Eq.(4-5) has the Schrodinger equation as a non-relativistic limit, and 
the Dirac equation as a one-body limit if one of (mi, 7712) — > 00. 

For the two- fermion state |2), Eq. ( |3-8| ) , to be simultaneously an eigenstate of J^, J3, and 
parity, the "bispinor" F = [Fij] must be of the form 



for 



Fir) 



Fir) 



_ 1 
r 

_ 1 

r 



isiir)tp^ir) + is2(r) V7°(f ), 
Mi(r) (/?~(f) -I- U2ir) (/3"'"(r), 

parity eigenstates, and 

■ isi(r)93-(f) + is2(r)(/3+(f), 
Ulir)^p^ir) +U2ir)(p^ir), 



tiir)(p (f) +t2ir)ip+ir) 
it!i(r) (f^ir) + if2(r) (f^ir) 



tiir)^^if)+t2ir) (^O(f) 
it)i(r) (f^if) + iv2ir) f'^ir) 



(6-1) 



(6-2) 



for —(—1)"^^^ = (—1)"' parity eigenstates. 

The 2x2 angular "bispinor harmonics" ip^ir), ip'^ir), ip^ir), ip~ir), for given total angular 
momentum quantum numbers J and rrij = M, are 



^^(f) = -^Yfir) 



-1 

1 



(6-3) 
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/(f) : 

and 



1 



^2J{J + l) 
1 



M-l 
J ' 



V(J-M + i)(j + M)y 
-My/, 

V(J + M-i)(j + M) 1 , 



- v/(J + M + i)(J-M) 



V(J + M)(j-M)y^^ 



J-1 

M+1 



V2J(2j- 1) L v/('^ + Af)(J-M)yf 1, vi-^- - i)(J-M)yj^+ 



(f) 



v/2(J + l)(2J + 3) 
v/(J - M + 1)(J - M + 2)yj^7 1 

J+1' 



(6-4) 
(6-5) 

(6-6) 



+ M + i)(j -M + i)y/ 



- v/(./ + M + 1)(J - M + i)yj^i 
V(J + M + i)(j + M + 2)yj*f;+^ 



We note that is antisymmetric and 99*^'^ are symmetric matrices. Furthermore ip^ , and 
correspond to opposite parity because Yjf {—r) = {—\)^Yjf{r) and f^", tp^ have L = J whereas 
y?^ have L = J it 1. These four bispinor harmonics form an orthonormal set, in the sense that 

/ df Tr((/7.^ ipj) = 5ij, where i,j = A,0,+, — and the integrations are taken over the entire sohd 
angle. 

The eight radial functions in the bispinors ( |6-lD and ( |6-2D are solutions of the coupled ra- 
dial equations that are obtained by substituting (|6-1|), (| 
coefficients of the four independent bispinor harmonics. 

We make use of the following identities in carrying out the radial reduction: 



into equation (4-5) and equating the 



<T'pf(r)u?(r) = — i— a ■ r ipir) H — f(r)cr • rcr ■ I (p(r), 
dr r 



(6-7) 



where (p{r) is any 2x2 bispinor harmonic, /(r) a radial function, and f = — , and I = r x p = 

r 

— ir X V. In addition, we note the following useful properties of the above bispinors harmonics: 



<T ■ r If 







where 



A 



I J + 1 
2 J + 1 



cr ■ r ip 

a ■ I 99 

a • I if'^ 



B 



Dip" + Aif'^ 

-(J + 2)99-, 
(j-i)(^+. 



and 



C = JJ{J + 1) 



(6-8) 
(6-9) 
(6-10) 
(6-11) 
(6-12) 
(6-13) 



(6-14) 



2J + 1 

It is evident from equations (6-1) and ( |6-2| ) that, in general, eight coupled radial equations are 
obtained, for arbitrary J > 0. 



7 Radial reduction of the two-body equations for J = states 



For the J = states, namely the 0^ (^'S'o) and O^(^Po) states, only two linearly independent bispinor 
harmonics arise, namely Lp^ and (p~ (equations ( |6-3| ) and (16-61) ), and so S2 = ^2 = 1*2 = 1^2 = 

(Here, as elsewhere, we give in brackets the nonrelativistic limit 



in equations (62I) and (^ 
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designation, "^^^^Lj, corresponding to the state in question.) Thus there is only one set of four 
coupled radial equations for each of 0^(^5o) and 0"^(^Po) states: 

K K 

(m+ + V{r) - E)s{r) - t'{r) - —t{r) - u{r) - — n(r) + ^V{r)v{r) = (7-1) 

K K 

(m_ + V{r) - E)t{r) + s'{r) s{r) + v'{r) v{r) + r]V{r)u{r) = (7-2) 

K K 

(-m„ + V{r) - E)u{r) + s'{r) s{r) + v' {r) v{r) + riV{r)t{r) = (7-3) 

K K 

{-m+ + V{r) - E)v{r) - t'{r) t{r) - u'{r) n(r) + £y{r)s{r) = (7-4) 

ds 

where m± = mi±m2, s' = —, the potential V(r) = —air [a = |gig2|/4vr), and E is the eigenenergy 

dr 

(two-particle bound-state mass) to be determined, while if = 1 (.^ = 2, r/ = 0) for the 0~(^S'o) states 



and A- = -1 (^ = 0,7? = 2) for the 0+(3Po) states. As shown in Ref. |], Eqs. (jT^) - (O) ^^ve 
the expected Schrodinger nonrelativistic limit, and the Dirac one-body limit. 

We note that the case with = r/ = corresponds to the simplified model without transverse- 
photon interactions, that is the 'Coulomb QED' model of Ref. [^. Similarly for = r/ = 0, if 



the potential is V{r') = , and the sign of the potential is reversed in (7-2) and (7-3) we 

VK r ^ — ^ ^ — ^ 



recover the 0^ radial equation of the Yukawa model discussed in Ref. ||3^, for which the inter 
fermion interaction is via a (massive or massless) scalar mediating field. 

We should point out that equations ( f7-l| ) - ( |7-4D , like the Dirac equations, have both positive 
and negative energy solutions. Indeed, in this two-body case, there are solutions of four types: 
E ~ mi + 7712, — —7771 + 7772, E ~ 7771 ~ "^-2 and E ~ — 777i — 7772, ^s cau be sceu most easily from 
the a = case. Of these, two are positive energy and two are negative energy solutions. 

Since we do not have analytic solutions for the eigenenergies of the present QED case, it is 
useful to illustrate this phenomenon on the scalar Yukawa (or "Wick-Cutkosky" ) model, in which 
scalar particles interact via a massive or massless mediating field. For such a scalar model, analytic 
expressions for the two-body bound state energy eigenvalues are available in the massless-exchange 
case [01: 



E 



\ 



where a is the effective dimensionless coupling constant, analogous to the fine structure constant 
of QED, and n is the principal quantum number. The it in Eq. (|7-50 correspond to two segments 
of a distorted semicircle. The upper branch of this distorted semicircle corresponds to the upper 
(positive) sign in ( [7-51 ). It begins from E = mi + m,2 at a = (indeed, E = mi + 7772 — |"^r(^)^ — 

^777^ (1 + (^)^ + • • • , for ^ <C 1), and decreases to a £^ = m^ + 7773 at the critical value of 
a = 77, beyond which E ceases to be real, and the wave functions cease to be normalizable. The 
lower branch, by contrast, begins from E = \mi — m2\ at a = and rises monotonically to meet 



the upper branch at the same critical point E{^) = ym^ + 7773 . These |777i — 7772I type bound 
state eigenenergies do not have the correct Balmer limit, since for this branch 

1 / 777x7772 \ / « \ ^ 1 / 7^17772 \ / 77717772 \ /o^^ 



E=\mi-m2\ + -i-. [- I ^^1-7 " +•••' M 

2 V |"7l — 7772 1 / \n / O \ 1 7771 " W-2| / V ("^1 — 1^2) / \n , 
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for mi / 771,2, but E = m ^(^) + |(^)^ + • • for mi = m2 = m. Thus, this "mixed energy" 
E ~ I mi — bound-state spectrum must be regarded as unphysical. There are also negative- 
energy solutions of the E ~ —mi — m2 and E ~ — |mi — m2\ type, but they are not bound-states, 
since the potential effectively reverses sign for the negative-energy solutions (as happens also in the 
Dirac-Coulomb case) . The same type of behaviour of the energy spectrum is observed in another 
analytically solvable case, namely a fermion and a scalar particle interacting via massless scalar 



quantum exchange |31|] Thus, we expect that the energy eigenvalue spectrum of (|7-1| ) - ( [7-4]) will 



be qualitatively similar to that of the scalar exchange models just described. 

We have not been able to determine solutions to the coupled radial equations (7^1)- (7^4) in 
terms of common analytic functions. It is of interest, therefore, to consider the properties and 
general behaviour of the solutions before commencing with numerical solutions. 

In analogy with the scalar model just described, and with the Coulomb QED case B we expect 



that, as a increases, the eigenenergy spectrum E{a) of equations ( |7-lD -(7-4) will have a qualitative 
behaviour similar to that of the Dirac spectrum, namely that E{a) decreases monotonically from 
E(a = 0) = mi + m2 until a hits a critical value Oc, beyond which E{a) ceases to be real. It is 
possible to infer the value of Uc by considering the ultra-relativistic limit, p — > oo, in which case 
we can neglect the masses mi and m2, and seek solutions of ( |7-lD -( [7^ ) with E = mi = m2 = 0. 
(This approach, when applied to the one-body Dirac-Coulomb case, yields the correct critical 
values ttc = \k\ = \j + 1/2|.) In this ultra-relativistic approximation, equations (7-l)-(7-4) have 



the solutions t = u, s = v, \t\ = \s\ = 1 (i.e. F oc -) with = 4i^^/(l + ^)(1 + rj), which gives 

r 

Oc = 2/-v/3 = 1.1547... for all 0^ states. Note however, that this result does not mean that the value 
of the two-fermion rest mass E at Oc is necessarily the same for the 0~ and 0^ states (certainly, 
such is not the case in the one-body limit). Note, also that the result, Oc = 2/\/3 for 0^ states 
is independent of the masses, that is we expect it to be the same for all finite mi/m2. The value 
Uc = 2/\/3 is different, and somewhat larger, than the known one-body limit (Dirac-Coulomb) 
value of = 1 for I'^l = 1 states. Also, this value is much smaller that the value Uc = 2, which is 
obtained for the two-fermion Coulomb QED case (where i = rj = {)) for 0^ states. 

Ol |'?l'?2| 

For the Coulomb potential V = , where a = — - — , it is often convenient to rescale the 

r 47r 

radial variable, that is to let p = r/a, where a is a suitable scale parameter. For example, the radial 
functions s, t, ti, v have the large r (negligible V and K/r) behaviour s ~ e"'', etc., for positive 
energy J = bound states, where a is given by 

1 [m\ - E'^][E'^ - mV\ 1 2 ^ 



AE^ or -2=m -[^^ if mi = m2. (7-7) 

Eq. (7-7) implies that a is positive only for \mi — m2\ < E < mi + 7712, which means that the 



bound state spectrum must lie in this domain (c.f. (|7-5D ). From this, and in analogy with the 
scalar model results, we can infer that the critical value E{ac = 2/\/3) lies between E = 777+ and 
E = |777„|, and likely closer to the former rather than the latter. 

With the rescaling p = r/a, equations ( |7-l| )-( [7^ ) become modified slightly, in that r is replaced 
by ap in all of them. For purposes of numerical integration of the radial equations the scale 
parameter a can be chosen to be anything that is convenient, be it that given in Eq. ( [7-7|) , or a = 1 

or a = — , or whatever. 
pa 

For a power series analysis of the radial equations it is useful to make the replacement s = se~P, 
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etc. Assuming solutions of the form 

s = p^lao + aip + a2P^ 

2 



t = p<[bQ + hip + h2p' + ■■■], 

U = p'^[co + Cip + C2P^ H ], 

V = p^[dQ + dip + d2p^ H ], 



(7-8) 

(7-9) 
(7-10) 
(7-11) 



we find , upon substitution into the radial equations for s, t, -u, v and equating coefficients of powers 
of p'^^^~^ , that the coefficients aj,bj,Cj,dj must satisfy the following recursion relations: 



^adi, = 0, 



(7 — K + v)ay — 5ay-i + a{m^ — E)b^^i — ab^ + — K -\- v)dy — 5d^^i — rjacy = 0, 

(7 — + v)ay — 5ajy_i + a{—m^ — E)cu-i — ac^ + {'y — K + i')dy — 5d^^i — rjoby = 

+ K + u)bu - Sbu-i + {j + K + v)cv - 6cu^i + a(m+ -|- E)dv~i + ady + ^aa^ = 0, 

where (5 = 1. If 5 = then ( |7-12| )-( [7-15D are the recursion relations for the power series represen- 
tations of the functions s(r), etc., rather than for s(r), etc. 

For = 0, and bearing in mind that a_i = 6_i = c_i = d^i = 0, (7-12)-( [7-15| ) yield four 



(7-12) 

(7-13) 
(7-14) 
(7-15) 



coupled homogeneous equations for the parameters oq, 60 j 
non-singular) solutions only if 



7 



1 



Co, do, which have non-trivial (and 



(7-16) 




for the J = states, for any values of mi, m2, whereupon 



do _ feo _ Co 



(1+Oa _ 2(7 -i^) 



Go 



ao ao 



2(7 + K) (l + 77)a 



(7-17) 



The condition (|7-16| ) implies that the radial equations have real bound state solutions of the 
form ( |7-8| )-( [7^TT| ) only for a < for any values of mi and m2- This, in turn, implies that etc < 

for the 0^ states for any (finite) values of mi and m2, in agreement with the ultra relativistic limit 
discussed above. This condition for bound states, a < is additional to the one that follows 



from Eq. ( |7-7| ), namely that |mi — m2\ < E < mi + 1112- 

The recursion relations ( [7-12D -( ^-15| ), with ( 7-16 ) and ( 7-17 ), can be used to generate the power 
series form of the solutions of Eqs. ( [7-8D -( 7^TTD . These series converge in the domain r a/m-^-, 
as discussed below and in Such a series can be used, for example, as a starting procedure for 
the numerical integration of the radial equations ( |7-lD-([7^. 

Unlike in the Dirac case, the recursion relations (|7-12| )- (|7-15| ) do not admit power series solutions 
of the form ( [7-8| )-( |7^TlD , which terminate at the same power, say = n' , so that a„'_|_i = 6n'+i = 
Cn'+i = d„/_|_i = 0. In particular, the ground state solution is not of the simple form 



aop'^, t = bop'^', u = cop"', 



dop^ 



(7-18) 



as it is for the two radial Dirac equations. This is perhaps to be expected, since in the Dirac case 

there are only two functions, say s and t, and four unknowns to be determined, namely — , 7, a 

ao 
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and E. Since the two coupled radial Dirac equations yield four equations (the coefficients of p"^ and 
of it is not surprising that a solution is obtained. In the present case, we have four coupled 

radial equations (7-1)- ([7-41) , which yield eight equations (the coefficients of p'^ and of p'^~^) to be 
satisfied by the six unknowns of the proposed solutions ( [r-lS| ), namely bo, cq, do, 7, a and E. Thus 



the system is overdetermined and no solution of the form ( 7-18| ) is possible. This situation persists 
for any solution of the form ( |7-8[ )-( [7^Tl ) where the polynomials all terminate at the same degree. 
Therefore, we shall solve the radial equations (7-l)-( [7^ ) numerically. 

Equations (T^l) - (|7-4| ) are not independent. Indeed, elementary manipulations of these equa- 
tions, namely subtracting Eq. ( |7-4p from Eq. ( |7-1| ) and similarly Eq. ( |7-3D from Eq. ( |7-2| ) show 
that 

E -1714^ - {I - ^)V{r) E -m^ - {l-ri)V{r) 



v[r 



-s{r), u(r) 



[I - ri)V{r) 



t{r). 



(7-19) 



E + m+- {l-i)V{r) ^' ' E + m^ 

Thus, the number of equations can be reduced from four to two. 

We introduce the auxiliary functions /(r) = s{r) + v{r) and g{r) = t{r) + u{r). Then, adding 
Eqs. ( fM]) to Eq. ( |7^ and Eq. (|^2|) to Eq. (j^3|), and using Eq. ( [7-19D yields the equations 



nr) = -f{r)+Wg{r)g{r), 



-9'[r) 



K 



9{r) + Wf{r)f{r) 



where 



Wf{T) = - 



E - V^{r) 
E-V^{r) 



(mi — 
E-Vr,{r) 

(mi -I- m2)^ 
E-V^{r) 



(7-20) 

(7-21) 
(7-22) 



and where 



V^ir) = {l + OV{r), V^{r) 

VJr) 



{l + r))V{r) 



il-OV{r), 
(1 -r/)y(r). 



(7-23) 



For the present QED case in the Coulomb gauge, for the 0~ states (for which ^ = 2,?] = 0), 
Vf) = Vri = y while = 2>V and = —V. In this case we see that Wf{r) is singular at 
ri = a/E = a'^(mr/E){l/mra), where m^ is the reduced mass and l/m^a is the reduced Bohr 
radius. This singular point is quite close to the origin (in units of the reduced Bohr radius) for 
small a. The appearance of this singularity may signal difficulties in the numerical determination 
of eigensolutions of the equations ( 7-201) by standard "shooting" methods. For the 0"'"(^Po) states 
(for which = 0,1] = 2), the singularity at ri = a/E occurs in Wg, but only if mi 7^ m2. Thus for 
the equal-mass O"*" states, equations ( |7-20D have only the usual 1/r singularities at the origin, and 
are amenable to solution by standard methods, as discussed below. 



8 Perturbative determination of the relativistic correction to the 
two-body eigenenergies for J = states 



Equations (|7-2C| ) can be written in the matrix form 



H\ 



where 



H 



e-Wf 
A — K 

dr r 



dr 



K 

r 
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and where e = E — (mi + m2). If Wf is replaced by W\ 



nr 
/ 



V and Wq by Wg'^ = 2fi, where 



H = mi?7i2/(mi + 1712), then Eq. ( |7-20D , or ( ^-1| ), is equivalent to the radial Schrodinger equation 
The first-order correction to the non relativistic energy enr = 



Ae 



— ^/xa^ is then given by 



nr |9n 



-2) 



If we expand the coefficients Wf and Wg (Eqs. ( [7-10| ) and ( 7-11 )) in powers of V/rrii, and keep 
only the lowest-order terms, we obtain 



+ W^ 



Wf 



Wg^ 



Wf 



1 



-(enr - V^f 



2m_(_ 



-3) 



(8-4) 



where 5 = This leads to the following O(a^) correction to the non-relativistic energy for 

the J = states: 



Ae 



1 



2mj 



2E„,(i-fl(v> + (i-a"{i'"> 



enr((l)) + ^((l+^) + (l-r?)^')((V)). 



(8-5) 



We use the notation {X) = (/nr|-'^|/nr)/(/nr|/nr) but {{X)) = {gnT\X\gnr)/{fnr\fm:)- 

For the 0~ (n ^^o) states (for which K = 1 and ^ = 2,r/ = 0) this formula gives for the lowest 
order relativistic correction the result 



Ae (0- (n 



1 /3 



n 



8 mi 



1 /i 
- + 2 — 

2 111+ 



(8-6) 



m2 



m case. This does not agree with the 



which becomes — ma^ in the equal mass mi 
64 ^ 
21 

known Positronium value of El]- This is not surprising, since the Breit equation, without 

modification, is known to give the incorrect fine structure for Hydrogen and positronium. Brown 



and Ravenhall argue that the reason for this, is the mixing of positive and negative energy 
one-particle states (which arises, in our formalism, because of our use of the "empty" vacuum 
(|3-4|)). This difficulty of the unmodified Breit equation is discussed in various works {e.g. Refs. 
p^ , |33| , [1^). The modification that is needed to bring the result into agreement with the observed 
fine structure of H or Ps is to subtract off the expectation value of the operator p3| , ^] 



H' 



a 



4m+r^ 



(3 — 2<Ti • Cr2 + 0"ir.(T2r 



air = o-i ■ r/r, 



(8-7) 



where, in this equation, we use the notation of |T^]. The expectation value of (8-7) (with respect 
to the non-relativistic eigenfunctions) is 



{H') 



n^{2J + l 



-fs = /^a"^ 



4/5 



n3(2J + l)m+' 
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where /s = 1 for the singlet S = states, while fs = 1/4 for the triplet S = 1 states with J > 0, 
but fs = for the triplet states with J = (see, also, section 10 below). The expression (8-8) gives 

the value —ma^ for the equal-mass ground state, which, when subtracted from the 'Breit' value of 
11 



-ma 



64 
subtract 



gives the expected positronium result 



21 



-ma 



More generally for arbitrary masses, if we 



from the expression (^-5|) , we obtain the "corrected" result 



8 m-4_ 



1 



1 

2 

2 m 



(8-9) 



This same result ( ^-9| ) was obtained previously for the n = 2 state by Darewych and Horbatsch 
[p^]b, who used a perturbative approximation on variationally derived equations. 

Somewhat surprisingly, the unequal-mass O(a^) corrections for arbitrary states seem to have 
been worked out fully only relatively recently. We refer to the work of Connell |p5[, who used 
a quasipotential formalism based on the work of Todorov [^ ], and of Hersbach [07]], who used a 
formalism based on a relativistic generalization of the Lippmann-Schwinger equation due to De 



Groot and Ruijgrok [|38[. Our corrected expression (18- 9|) agrees with the results of these authors. 
(The O(a^) corrections for Hydrogen and muonium quoted in standard references are expansions 
in mi/ma (e.g. ||], @.)) 

For the nO"*" (^Pq) states (for which K 



-1 and ^ = 0,?7 = 2) equation (8-5) gives 



Ae (n 0^ 



'Pa 



^a 



8 m_|_ / 



+ 



2 

3 m^ 



S-10) 



which does agree, in the equal mass case, with the Ps values for all the n O^('^Po) states, as well 
as with the unequal- mass expressions of Connell [^] and Hersbach |37| for these states. This 
agreement implies that the "correction" {H')^^ vanishes for the ^Pq states, as indeed it does. 

We might note, in passing, that formula ( ^-5| ) gives the correct O(a^) results for the Coulomb 
QED (C = = 0) case for which the W coefficients are non-singular for r > 0. 



9 Numerical solutions for some J = 0+ states 



In the case of equal masses mi = m2 = m the radial equations (7-20) for J = 0^( Pq) states are 
free of singularities. Thus the boundary value problem is well posed, and it can be solved by means 
of a standard numerical ODE-solving procedure. We solved it by the "shooting" method using the 
Maple Runge-Kutta programme. 

The corresponding perturbative spectrum (c.f. ( 8-lOD ), 



a 



2 



Elm = 2-—. + a^{ —— ^ ^ , (9-1) 

' 4n2 \ 64n4 Sn^ J ' ^ ' 

agrees with the orthopositronium spectrum since the contribution of the extra terms (c.f. ( |8-7D ) 
(caused by positive-negative energy mixing) vanishes in this case. 

In the Table 1 the numeric and perturbative results are presented for the lowest-energy 0+ state 
(i.e., J = 0, £ = 1, n = 2) for different values of a < = 2/^3 ^ 1.1547005383792515. 

Note that the perturbative O(a^) results are unreliable for a > 0.5. Our numerical results 
suggest that, for this n = 2, 0"^ state, E{ac = 2/\/^) ~ 1.29974 m, which is smaller than the scalar 
theory value E{ac = 2)/m = V2 = 1.4142.. (see below (|7^)). 
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Table 1: Values of E/m for the n = 2, O^(^Po) state (mi = m2 = m). 



(X 


Perturbative 




Numeric 




(equation 




'''' 


(equations (7-20|)) 


1 /I "Vl 


1 
1 


QQQ 


QQR 


669 

\J\JiJ 


Q^'^9 


1 QQQ 


QQR fifiQ Q^'^9 


U.Ul 


1 
1 


QQQ 


QQ"? 


74Q 


fiQDS 

uyuo 


1 QQQ 


QQ"? 74Q fiQDS 


u.uo 


1 
1 


QQQ 




00\J 


799n 


1 QQQ 




n 1 


1 


999 


371 


Q07 


OOZl 


1.999 


371 886 




1 


982 


442 


220 




1.982 


028 02 


0.7 


1 


961 


950 


032 


5521 


1.957 


997 74 


0.9 


1 


929 


085 


449 


2188 


1.902 


4531 


1.0 


1 


906 


575 


520 


8333 


1.838 


781 05 


1.1 


1 


879 


098 


470 


0521 


1.688 


2317 


1.15 


1 


863 


256 


642 


6595 


1.436 


9434 


1.154 


1 


861 


924 


191 


0995 


1.355 


170 76 


1.1547 


1 


861 


689 


995 


0549 


1.301 


3199 


1.1547005383792 


1 


861 


689 


814 


8148 


1.299 


74 



There are, as we explained in section 7, > "mixed energy" solutions of the form E/m = 
a + O(a^), which are unphysical, because they do not have the Balmer non-relativistic limit. We 
do not list such solutions here, though they can be calculated readily enough in the same way as 
those in the Table 1. This unphysical branch of the n = 2 0^ state rises uniformly from zero at 
a = to join the physical branch of Table 1 smoothly at E{ac)- As mentioned previously, the two 



branches together resemble a distorted semicircle (c.f. Eq. (7-5)). 

Analogous results for the n = 3 and n = 4 0"'"(^Po) equal-mass two-fermion energies are 
given in Table 2. The qualitative behaviour of E{a)/m for these states is similar to that for 
the lowest such state [n = 2), except that the critical value of a increases with n, as it does 
in the case of the analytically solvable scalar model of Eq. ( [7-5| ). However, here we obtain 
ttc/n = 0.64987, 0.6024745, 0.4871785 for the n = 2, 3 and 4 states respectively, in contrast to 
the scalar model values adn = 1 for all n. 

Note that the critical value of the two-body mass, E{ac)/'m, increases with n, in contrast to 
the scalar model, for which E{ac)/m = \pl for all n. 

Figure 1 is a plot of the unnormalized reduced radial wave functions s(r), t(r) = ti(r) and 
v(r) (see Eqns. ( |7-lD -( [7^ )) in the case of equal massless, ra\ = m-2 = m, for the lowest-energy 
n = 2 0"'"(^Po) states, when a = 1. These wave functions are qualitatively similar to those obtained 
for these states in the Coulomb QED case [^. The "large" component s(r) is nodeless while the 
"small" components t[r) = u{r) and the "doubly small" one v{r) have one node. The node at the 
origin, r = 0, is a consequence of our use of reduced radial wave functions s(r), etc., rather than 
the actual s(r)/r, etc. Indeed, the wave functions behave at small r as follows: 

s{r) ~ v{r) ~ const • ar"' , t{r) = u[r) ~ const • (2/3)(7 -|- l)r'^, (9-2) 

where < 7 < 1 (see Eq. ( [7-1^ )). Thus the matrix wave function F{r) is singular, F{r) ~ r'^^^ , 
as happens also in the one-body Dirac equation with a Coulomb potential. Nevertheless, F{r) 
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Table 2: Values of E/m for the n = 3,4 O^(^Po) states (rn-i = m2 = m). 
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is normalizable for all a up to and including a = ac, at which point t(0) = n(0) = (2/3),s(0) = 
(2/3)w(0) 7^ 0, as can be seen in Figures 2 and 3. 

Figure 3 represents the excited n = 3 O^('^Po) state for the critical coupling strength ac = 2/\/3. 
In this case, the wavefunction s(r) has one node while t{r) and v{r) have two nodes. This behaviour 
differs from that found in CQED |^] for the same case (where the number of nodes was two, one 
and three, respectively). In the n = 4 O^(^Po) case (Figure 4) each of the functions s{r), t{r) and 
v{r) gets one more node . This tendency likely continues for higher values of the quantum number 
n. 
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Figure 1: Reduced radial wavefunctions for the lowest m\ = = m n = 2 0+(3Po) state for 
a = 1, E/m = 1.838781. s{p): full curve; t{p) = u{p): broken curve; v{p): chain curve, p = r/a, 
where a = 2.542291(l/m). 
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1 2 3 4 5 6 r/a 

Figure 3: Same as Figure 1 but for the excited n = 3 0"'"(^Po) state, with a = etc = 2/\/3, 
E/m = 1.870423 and a = 2.824148(l/m). 
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Figure 4: Same as Figure 1 but for the excited n = 4 0+(2Po) state, with a = 1, E/m = 1.974982 
and a = 6.342094(l/m). 
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10 Radial reduction for J > states 



For states with J > 0, the eigenstate problem reduces to a set of eight first-order differential 
equations for the functions si{r) . . . V2{r) and the energy E {c.f. Eqns. (6-1) and ( |6-2| )). It is 
convenient to present this set in the following matrix form. Let us introduce the 8-dimensional 
vector-function: 

si(r) 



X(r) 



S2W 



(10-1) 



Then the set of radial equations reads: 

nX{r) = |h-^ + U(r)| X(r) = EX{r), 

where the 8x8 matrix U(r) has the following structure: 

U(r) = M + [G-a(l + S)]/r. 
Here I is the unit matrix, M is diagonal, 



(10-2) 
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(10-4) 



m± = mi lb 777-2, and the form of 8x8 matrices H, G, and S depends on the parity P: 
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(10-7) 
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for P = (-) 
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(10-9) 



(10-10) 



for P = (—)"', where A, B and C are defined in Eq. (|6-14D . Due to the properties 

H^ = -H, U^ = U, 
the "radial" Hamiltonian Ti \s a, Hermitian operator with respect to the inner product 



oo 

(Y|X)8 = j drY\r)X{r), 



:io-ii) 



(10-12) 



where the subscript "8" denotes the dimensions of the vector-functions X, Y. 

In the subsequent reduction of the set ( 10-21 ) one can use the fact that rank (H) = 4 (for either 
parity). Thus one can reduce the number of differential equations from 8 to 4. We perform this 
reduction in a way that ensures, as far as possible, the Hamiltonian structure of the equations. 

First of all we perform the orthogonal transformation: 



X(r) = EX(r), H = EHE 
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(10-13) 
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(10-14) 
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(10-15) 



It preserves the inner product ( 10-12| ) and reduces the equations ( 10-2| ) to the form: 
7YX(r) = (h-^ + U(r)l X(r) = SX(r), U = EUE^^ 



dr 

It is convenient at this stage to express the 8-dimensional vectors and matrices 
4-dimensional blocks: 



Xi 

X2 



V(r) = y]{r)-E\ 



Vii(r) Vi2(r) 
V2i(r) \/22(r) 



etc. 



Then the matrix H takes the form 

Hii 
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11 



JO 




" -1 " 


J 


J = 


1 



and the set of Eqs. ( |10-16 ) becomes 



HnXi(r)+Vn(r)Xi(r)+Vi2(r)X2(r) = 0, 
V2i(r)Xi(r)+V22(r)X2(r) = 0. 

The set ( 10-20| ) is purely algebraic. It permits us to express X2(?') in the terms of 

Mr) = -%2{r)^2i{r)Mr)- 
Substitution of ( |10-21| ) into ( |10-19 ) yields the closed set of four first-order differential 



C{E)Mr) ^ 



iHn|: + \A/(S)}xi(r) = 0, 



(10-16) 
in terms of 
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equations: 

(10-22) 



where 
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(10-23) 
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Here the upper sign corresponds to P = {—)'^~^^, and the lower sign corresponds to P = (— )'^. 

The operator C{E) in Eq. ( 10-22| ) is formally Hermitian, i.e., given E, it is Hermitian with 
respect to the inner product (• . . | • • .)4. But any two solutions of ( |10-22|) , Xi and Yi, corresponding 
to different values of the energy, E and E' , are not orthogonal. This is due to the nonlinear 
dependence of C{E) on E. Orthogonality can be instated by using the following definition of inner 
product: 

C{E') - C{E) 



((Yi|Xi))4= Yi 



(10-24) 



E' - E 

This inner product follows directly by substitution of Eq. ( 10-21 ) into Eq. ( 10-12| ). 

The set of first order equations ( 10-22| ) can also be expressed as a second order equation. For 
this purposeit is convenient to permute the elements of Xi by means of the matrix L, where 
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Then, in terms of the two-dimensional blocks 

Xi = LXi 

the equations take the form 
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(10-27) 
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where W 
*i: 



L W L ^ . Elimination of ^'2 leads to the following equation for the 2x1 vector- function 
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(10-30) 



Eigenstates ^'i, corresponding to different values of the energy are orthogonal with respect to 
the inner product: 

T{E') - T{E) 



((^il^i)), 



$1 



E' -E 



^1 



(10-31) 



which also follows from the reduction procedure. 



11 Perturbative solutions for J > states 

The form of equations ( |10-3C| ) is convenient for examining the energy spectrum of J > bound 
states perturbatively to O(a^). For this purpose we introduce the dimensionless quantities, 

E — rrij^ m_ 
p = IJ,ar, A = 5 — , and o = , (H-l) 
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where fi = mim2/m^ is the reduced mass. We now perform a perturbative expansion in a of 
equation ( 10-3C| ). To order a"^ equation ( |10-30D takes the Hamiltonian form (hereafter we omit the 
subscript "1" of the wave function 



where 



V'i(p) 

^2{P) 



(11-2) 



(11-3) 



is a two-component wave function, 7i is the Hamiltonian divided by pa"^ and expressed in terms of 
dimensionless quantities ( 11-1| ), 



7^(0) 



d_ 

2 Id^ 



1 
P 



is the unperturbed (i.e., Oth-order) Hamiltonian, and 



|^p,A)|- + A.(p,A) 



(11-4) 



(11-5) 



is the perturbative correction to ( 11-4 ). The form of the symmetric 2x2 matrices J^, IC{p,X), 
M{p, A) depends on the parity. In order to obtain the energy spectrum to O(a^), it is sufficient to 
calculate the eigenvalues of A to O(a^), i.e., A ~ A^'^^ +a^A*-^\ where A^'^^ will be calculated exactly, 
while for A^^-* first order perturbation theory in is sufficient. Hence, the dependence of '^^^•'(A) 
on A is not crucial: to the accuracy required, it can be replaced byA'-'^^. In addition, the kernel of 
the inner product ( 10-31]) can be set to unity. 

In the case P = (— )'^^^ we have J = C^\, so that TC^'^^ is the (dimensionless) radial Coulomb 
Hamiltonian Hj with the angular momentum i = J, repeated twice: 
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The matrices IC{p, A) and M.{p, A) are 
1 
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The eigenvalues of the Oth-order Hamiltonian ( p.l-6| ), namely 



A(o) = -l/(2n2), 



n = 1,2,..., 



are two-fold degenerate, each with the two eigenstates 
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(11-9) 
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where \n, J) is a solution of the Coulomb problem Hj\n, J) = X^^^\n, J). Thus, the correction A^^-* 
must be calculated appropriately for the degenerate situation: 



^(1) 

Xl,2) 



All + A22 ± 



(All - A22)2 + 4A? 



12 



fll-lll 



where the matrix A is defined as follows: 



A 
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The mass spectrum 



2C(2J+l)n'l 
_ _(3+£f)C^+2 
32n4 4C^(2J + l)ni' 



(11-12) 



^11-13) 



obtained with the use of ( |11-9D , ( |11-11 ) and ( 11-12| ) coincides neither with the muonium spectrum 
found in |35, 37| nor (if rn-i = 1712) with the spectrum of parapositronium (see |18] and refs. therein). 
The reason, like for the J = states, lies in the use of the empty vacuum (3^4), which leads to 
single-particle states of positive and negative energy, and subsequently to the Breit equation with 



its spurious term H' ( |8-7] ) in the Hamiltonian |P3| , 34, 18 1, as discussed in connection with the J = 
states. We are going to show that the elimination of the contribution of H' from the spectrum 
leads to the correct result. 

First of all we transform H' into a radial representation. For this purpose we note that ipi = 
si + ^'i (~ si in the nonrelativistic limit), i.e., ipi contains only those components of F{r) ( |6-1D 
which are coefficients of "bispinor harmonics" 0^(f). Similarly, ■02 = —S2 + V2 (~ — S2 in the 
nonrelativistic limit), i.e., t/ji contains only coefficients of the "bispinor harmonics" (lP{f). Thus, 
for P = {—Y^^ parity states the spurious term (divided by /ia^) takes the following radial form: 
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where i, j = A, 0, and the corresponding matrix elements are 
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If we now use A = A — A', instead of A, in ( [11-llD and ( 11-13 ), we obtain the spectrum 
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^ ^a^/ ll + j2 2 J + 1 VI + 4C2(^2 I 
2^^4^j~8^5^ (2J + 1)C2 j 



(11-16) 



which coincides with the results of Connell []35|| and Hersbach for the parity (— 1)''^^ states. 
Thus, correcting for the spurious terms in the Breit Hamiltonian, we obtain the expected O(a^) 
results. 

In the P = (— I)'' case, the matrix J is not diagonal: 
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(11-17) 
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It can be diagonalized by means of the orthogonal transformation, using the matrix: 

n 

Then 
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so that the Oth-order Hamiltonian becomes: 
It possesses the doubly degenerate eigenvalues ( 11-9| ) with eigenstates: 
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The first-order correction TY^^^ ( [Ll-5| ), with the matrices 
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generates diagonal matrix elements only: 
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The energy spectrum ( 11-13 ) with a|J| = An, A|2j = A22 contains the contribution of the spurious 



(11-25) 



term ( |8-7] ). Again, we present this term in radial form: 
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where, for the present P = (— ) parity case, i,j = —,+■ The corresponding matrix elements 
are: 
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The substitution of A^^^j = An — A^ and A^2) = ^22 — ^'22 i^^to ( 11-13 ) yields the spectrum: 
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which coincides with the results of Connell [Esl and Hersbach for the parity (—1)"^ states. 
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12 Concluding Remarks 



We have studied a reformulation of QED, in which the coupled Dirac-Maxwell field equations are 
partially decoupled by expressing the mediating photon field in terms of the Dirac particle field, 
using covariant Green's functions. This allows us to reformulate the Hamiltonian of the theory 
so that the photon propagator appears directly in a quartic, nonlocal interaction term. We then 
consider a truncated model, in which there are no free (physical) photons. For such a model, each 
iV-particle segment of the Fock space of the quantized, equal-time Hamiltonian is an invariant 
space, that is, there is no coupling among the various A^-fermion segments. This is achieved by 
introducing an unconventional "empty" vacuum state. As a consequence, there exist exact few- 
particle eigenstates of the truncated Hamiltonian, which lead to Dirac-like two- and three-fermion 
wave equations. We show, in particular, that the two-fermion wave equation, in the Coulomb 
gauge, is just the Breit equation. 

For specific states, the Breit equation reduces to radial form, and then to Dirac-like equations 
for J = states, and to a coupled pair of Schrodinger-like equations for J > states. Perturbative 
solution of these equations yields a^-corrections to the nonrelativistic Rydberg spectrum, which 
do not reproduce the muonium spectrum as calculated by Connell fs^ and Hersbach [^] (nor 
the positronium spectrum in the mi = m2 case). The apparent reason for this disagreement is 
the mixing of positive- and negative-energy states, which is characteristic of the Breit equation 
[p^, 18 1 . However, agreement is achieved if we subtract the contribution of spurious operator 
(FtI ), which appears in the Breit equation (c./. Eqs. (|8^ ), (|1 1-161) , (|11-27| )). 



We have not been able to obtain analytic solutions of the radial equations. These radial wave 
equations have, in general, a singular point at ri ~ a/E' > 0, where E is the two-fermion bound 
state energy (rest mass) . The existence of such an "interior" singularity makes it difficult to obtain 
numerical solutions of the radial boundary value problem by standard methods. The only exception 
is the case J = 0^, mi = m2, for which the radial equations are regular, and which we have therefore 
studied numerically. (It is noteworthy that this is the case where the contribution of the operator 
(13) is zero.) 

Our numerical results for the equal-mass J = 0+ states, show that the dependence of the energy 
E on the coupling constant a is qualitatively similar to that obtained earlier for the Coulomb 
QED model (for which transverse-photon interactions are ignored) Q. For low a, the numerically 
obtained eigenenergies are in agreement with the result derived perturbatively. Thereafter, E{a) 
decreases monotonically to E{ac) > as a approaches a critical value Oc- We find that Oc = 2/\/3 
(in contrast to CQED value of Uc = 2). 
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